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ABSTRACT 

The breaking stress (the maximum of the stress-strain curve) of neutron star crust is 
important for neutron star physics including pulsar glitches, emission of gravitational 
waves from static mountains, and flares from star quakes. We perform many molecular 
dynamic simulations of the breaking stress at different coupling parameters (inverse 
temperatures) and strain rates. We describe our results with the Zhurkov model of 
strength. We apply this model to estimate the breaking stress for timescales ~ 1 s 
- 1 year, which are most important for applications, but much longer than can be 
directly simulated. At these timescales the breaking stress depends strongly on the 
temperature. For coupling parameter V < 200 matter breaks at very small stress, if it 
is applied for a few years. This viscoelastic creep can limit the lifetime of mountains on 
neutron stars. We also suggest an alternative model of timescale-independent breaking 
stress, which can be used to estimate an upper limit on the breaking stress. 

Key words: equation of state - stars: neutron - stars: interiors. 



1 INTRODUCTION 

Ions in the neutron star crust can form a Coulomb crystal 
which determines the crust's elastic properties. For example, 
the breaking stress <7b is the maximum stress, as a function 
of strain, that the crust can support. If larger stress is ap- 
plied to matter it will not remain in a static configuration. 
Neutron star crust is und er high pressure. As a res ult, voids 
or fractures do not occur (|Horowitz fe KadaufeOOS^ ) and this 
simplifies how the crust breaks. The crust will break if the 
local stress is larger then the breaking stress at at least one 
point in the crust. In this paper we mostly refer to break- 
ing stress (force per unit area). The breaking strain (corre- 
sponding value of the fractional deformation) can be esti- 
mated from the linear stress-strain relation, where the slope 
is given by the well known elastic constant or shear modulus 
jStrohmaver et al.lll99ll : iHorowitz fc Hughtolliool ) . 

There are a lot of models which associate breaking of 
neutron star crust with observational phenomena. First of 
all, crust breaking c auses pulsar glitches in the starquake 
model of iRudermanl ll 19691, Il99lll and in some recent mod- 
els (for brief review see IChamel fe Haensell l200sh . Second, 
some models of magnetar gia nt flares involve crust breaking 
^Thompson fc Duncan! l200ll ). Here the crust may need to 
be very strong to control the release of large magnetic ener- 
gies responsible for extremely energetic gamma ray bursts. 
Indeed, our results, see below, predict that the crust is 
very strong. Finally, the breaking stress limits the max- 
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imum size of "mountains" on neutron stars. Mountains, 
on rapi dly rotating stars, can efficiently em i t gravitationa l 
waves (|Ushomirskv. Cutler fc Bildstenl I2OO0I : iHaskelll l2007h 
that could be detected by pr esent large scale interferome- 
ters (|Abbott et al.ll2007l . |2008) . Gravitational wave emission 
could be especially important for low mass neutr on stars 
that can have large deformations |Horowitzl l2009h . More- 
over the balance of angular momentum gained from accre- 
tion, and radiated in gravitational waves, can contro l the 
spin period of some accreting stars (|Watts et al.ll2008h . 

For most terrestrial materials the durability (time be- 
fore breaking) r at a given stress is known to depe nd on the 
temperature and, of course, the applied stres s (e.g., Zhurkovl 
1 19571; iRegeb, Slutsker fc Tomashevskiil 1 19721 : [siutsked 120051 
ISlutsker et al.l 12007). In other words, the breaking stress is 
not just a constant, defined by the matter parameters (den- 
sity, temperature and composition), but depends on the du- 
ration of the stress - the matter can break at lower stress, 
if one waits long enough. In fact, the dependence of o~h(r) 
on r is logarithmic and we will refer to r as the timescale 
of the process. The aim of the present paper is to determine 
how Ob, depends on r and on temperature T for neutron 
star crust. We suggest a simple expression [see Eq. © pa- 
rameterized by Eq. (O below] for the durability of neutron 
star crust material at a given stress. This expression gives 
the lifetime of elastic deformations (mountains, for exam- 
ple) and one can easily extract the breaking stress for the 
timescale of the process of interest (a few years for pulsar 
glitches, for example). 

We also note the possible application of our results to 
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the physics of dusty plasmas. Such a plasma can form a 
Yukawa crystal, just like the ions in neutron star crust. 
Its breaking stress can be measured in the laboratory as 
a threshold stress which is needed to obtain flow. The pres- 
ence of such a t hreshold has be e n rece ntly experimentally 

demonstrated bv lGavrikov et al l <|2009h . 

Until recently bre aking stress studies (|Smoluchowskil 
Il97d ; lRudermanlll99ll ) have been based on analogies with 
terrestrial materials instead of accurate calculations, and the 
uncertainties were large. Recently, large scale molecular dy- 
namics (MP) simulations of the breaking stress were per- 
formed IIHorowitz fc Kadaull2009h . The breaking strain was 
found to be large ~ 0.1 and only slightly affected by impu- 
rities, defects, and poly crystalline structure (grain bound- 
aries). In the present paper, we extend these results with 
extensive MD simulations of breaking stress at different tem- 
peratures and strain rates. We use these results to extract 
par ameters of the Z hurkov mode l of st rength [see Sec. 12.21 
and lZhurkovl l| 19571 ) ; [Regel' et all (|l972T ) for details] and ap- 
ply this model to estimate the breaking stress for timescales 
of 1 s - 1 year, which are of most interest for applications, 
but much longer than possible for direct MD simulations. 



2 FORMALISM 

The neutron star outer crust consist of ions and elec- 
tro ns. The inner crust also contains free neutrons (s ee 



^ Haensel. Potekhin fc Yakovlevll2007l ;l Chamel fc Haensel 
2008). The electrons are degenerate and form a slightly po- 
larizing background which screens the Coulomb interaction 
between ions. We describe electron screening by a Thomas- 
Fermi screening length 



1/2 



(37r 2 n e 



-1/3 



(1) 



where n c is the electron number density, e is elementary 
charge and c is the speed of light. For simplicity, we assume 
the electrons are ultrarelativistic. 

In this paper we discuss a one component plasma, where 
all ions have the same charge number Z. The ions interact 
via a Yukawa potential 



cxp (-r/A c ) , 



(2) 



where r is the interion distance. The state of the ion system 
can be characterized by the classical coupling parameter 



aT 



(3) 



Here a = [3/(47rrii)] 1 ^ 3 is the ion sphere radius and T is 
the temperature in energy units. The ion number density is 
7ii = n c /Z because the system is neutral. 

In the ultrarelativistic limit, the ratio of the electron 
screening length A c to the ion sphere radius a depends only 
on Z: A c /a w 5.41/Z 1/3 . We use Z = 29.4 in all of our cal- 
culations. This is th e mean charge of the ions in the crust 
composition used by I Horowitz. Berry fc Brown! (120071 ) and 
iHorowitz fc Kadaul |2009l ). It corresponds to A e /a ~ 1.75. 
For such screening length the crystal is thermodyna mically 
stable at T > T m « 176.1±0.7 jHorowitz et alj|2007h . which 
is close to the melting point r m = 175 expected for the 



one component plasma in th e absence of electron screening 
jPotekhin fc Chabrierl I2OO0I ) . We suppose, that the break- 
ing stress depends only slightly on screening length and our 
results can be used for matter of any composition. 

A typical dynamical frequency is the ion plasma fre- 
quency w p , 



47rZ 2 e 2 rii 



(4) 



where mi = Am u is the ion mass. The ratio of plasma tem- 
perature T p = hw p to the temperature T characterizes the 
importance of quantum effects on ion motion. In our clas- 
sical MD simulations we neglect quantum effects (assuming 
h — > 0), but they can be important at the conditions of 
realistic neutron star crust and we discuss them in Sec. 13.21 



2.1 Molecular dynamics simulations 

To calculate the breaking stress we developed a parallel ver - 
sion of the YUKAWAMD code l|Horowitz fc Kadaul |2009| ) 
where the ion system is strained by deforming periodic 
boun daries. We use the velocity Verlet algorithm l|Verletl 
Il967h with a time step ~ 0.1/w p . The simulations are per- 
formed as follows: the original system is taken to be a per- 
fect body centered cubic crystal contained in a cubic box, 
that is aligned with the lattice planes. After thermalization, 
we start a shearing deformation of the periodic box bound- 
aries according to x — > x + ey/2, y — > y + ex/2, and z — > 
z/(l-e 2 /4). This deformation conserve volume of the simu- 
lation box. Here the strain e — vt increases with time at con- 
stant strain velocity v. Since the deformation is adiabatic , 
we define stress as a = d£/de (|Landau fc Lifshitd Il986h , 
where £ is internal energy at unit volume. The derivative 
is calculated numerically as <r(e) = [£ (e) — £ (e — Ae)] /Ae 
with Ae = 0.002. The breaking stress <Tb is the maximum 
stress obtained during the shear simulation. The correspond- 
ing strain we refer to as the breaking strain eb- 

To decrease the simulation time we prepare a well ther- 
malized system at strain e = 0.05 and use it as an initial con- 
figuration for deformations with different v. For the longest 
runs (with lowest v) we start from a configuration at larger 
strain. The initial parameters of this configuration are ex- 
tracted from other simulations with larger v. We also use a 
cutoff radius r cut to omit interactions at r > r cut and in- 
crease computation speed. Most of the simulations are done 
for 9826 ions in a periodic box and r cut = 13A C . This large 
r C ut valu e is needed to accurately r eproduce the shear mod- 
ulus (see lHorowitz fc Hu ghto 2008). We discuss these simu- 
lation parameters in more detail in Sec. [3] 



2.2 Zhurkov model of strength 

The Zhurkov mode l (see e.g. IZhurkovl 1 19571 ; iRegel' et all 
Il972l; ISlutskeij|2005l ) is based on the kinetic conception of 
strength. The main assumption is the following. The break- 
ing event occurs due to thermal fluctuations. To achieve 
breaking, a fluctuation should have energy U. This energy 
threshold is reduced by the product of the stress a times 
an activation volume V. The probability of a breaking fluc- 
tuation is oc exp [—(U — aV)/T]. Following this model, the 
matter will break at stress ah, if it is applied over the time 
interval r 
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Figure 1. Dimensionless breaking stress <5b versus dimensionless 
shear speed v/u) p . The symbols are MD data (dots, circles, filled 
triangles, triangles, filled rhombus, rhombus, filled squares and 
squares are for T = 180, 200, 250, 300, 400, 800, 1600 and 3200 
respectively). The lines are predictions of the Zhurkov model, Eq. 
iffjl. with parameters given by Eq. (gj. 



to exp 



U-abV 
T 



(5) 



where to is a typical timescale for ion motion, that we take 



to be ro 



Denning dimensionless threshold energy 



U = Ua/(Z 2 e 2 ) and stress a = a/(mZ 2 e 2 /a), Eq. © can 
be rewritten in the form 

r = — exp (UT - a b NT) , (6) 

where N — Vrii is the number of ions in the activation vol- 
ume. We use U and N as fit parameters to reproduce MD 
data. Equation ([6]) is written for static stress, but our MD 
simulations are dynamical. The strain is linearly increasing 
with time, this produces a stress that also increases, almost 
linearly, with time a « fie = fivt. Of course, the linear de- 
pendence of a on t is not correct after the breaking event. 
Here fi w a b /e b is the corresponding elastic constant. Equa- 
tion Q can be easily gene ralized for such deformation (see 
ISlutsker fc Aidarov I (| 19831 ). for the example). The strain ve- 
locity v and breaking stress are related by the equation 

- = T^exp(-L7r + a b 7Vr). (7) 
cop NFa b 

The stress is lower than a b during most of the deformation 
time. As a result, the time before breaking r = e b /v is larger 
than t given by Eq. © by a factor of NFa b . 



3 RESULTS 

Our MD simulation results for the breaking stress of a sys- 
tem of 9826 ions are shown by symbols in Fig. [T] The data 
are fitted by Eq. (0 with parameters U and N(T), 

500 



This fit is shown as lines in Fig. \T\ The e b in Eq. was 
taken from spline interpolation of MD data over v/ui p . In 
fact, a b depends on e b only logarithmical and one can set 
£b = £b(r) (and neglect the dependence on v/ujp) in Eq. 
to reproduce the plot. There is good agreement between 
the fit and MD results, except for v/uj p >2x 10 -5 . We sup- 
pose, that here the deformation velocity is too rapid and the 
breaking stress increases because there is not enough time 
for the system to break and relax the lattice. We exclude 
these points from the fit. 

If Eqs. (|6]) and JH} are formally applied for F below 
the melting value F m , we find that the crystal has a non- 
vanishing durability at fixed stress until r > 150. T h is is i n 
qualitative agreement with the results of Daligault (2006), 
who found Coulomb crystals to be metastable for T > 150. 

3.1 Dependence on system size and cutoff radius 

To study finite size effects we have performed a simulation 
with 250000 ions at T = 800. For a strain speed v/uj p = 
6.25 x 10 -6 we obtain a breaking stress a b = 0.0187. This 
is only 4% larger than the value 0.0180 obtained for a 9826 
ion simulation. We conclude that finite size effects are not 
very large. 

To study the dependence on cutoff radius r cut we per- 
form a set of additional runs with r cut = 10A o and r cut = 
16A e . The value r cut = 10 A c is too small (for F = 800 
the breaking stresses are systematically 10% larger than for 
r C ut = 13A C ), but results for r cut = 13A e and r cut = 16A e are 
almost the same (within a few percent statistical accuracy). 
We conclude that r cut = 13A C is large enough. It is easy to 
understand why so large a value of r cut is needed. Let us 
estimate the interaction energy U C ut of an ion with other 
ions at distances r > r cut . Assuming the ions are uniformly 
distributed for r > r cut one obtains 



Ucut — 3 



Z 2 e 2 A c (r cut + A c ) 



exp(— rcut/Ac). 



(9) 



U = 0.366, N 



149 



+ 18.5. 



(8) 



For r cut = 10A e and A c = 1.75a we obtain J7 cut w Z 2 e 2 /220a 
which is small compared with the Coulomb energy, but for 
T > r m w 176 it is of the same order of magnitude as 
the thermal energy, which is crucial for the breaking of the 
crystal. For r cut = 13A e the energy f7 C ut ~ Z 2 e 2 /3400a, 
which is small compared with thermal energy up to T < 
1600. 



3.2 Correction for quantum effects 

In our classical MD simulations we can not properly in- 
clude quantum effects, such as zeropoint vibrations. But 
such effects can be important for applications, since the crust 
temperature can be less than the plasma temperature. At 
T < T p the breaking event takes place because of subbar- 
rier tunneling, but not overbarrier thermal fluctua tions as in 
a classical crystal (see ISlutsker fc Aidarov 1 11983. for exam- 
ple). To include quantum effects we make a simple assump- 
tion that the breaking stress mainly depends on the root- 
mean-square displacement of the ions. A similar idea was 
susg ested by ISalganikl l|l97fj) to describe polymers break- 
ing and the results were shown to be in a good agreement 
with experiments with boron samples ( Slutsker fc Aidarov I 
Il983h . Let us introduce a renormalized coupling parameter 
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Figure 2. The ratio of the breaking stress <r b to the electron 
pressure P c 4.4 X 10 26 dyncm - 2 (left axis) versus timescale of 
the stress r at five temperatures T = 1 X 10 6 , 2 X 10 7 , 5 X 10 7 , 
1 X 10 s and 2 X 10 8 K for 56 Fe matter at the density 10 9 gem -3 . 
The solid lines are predictions of the Zhurkov model, Eq. 10, 
with parameters given by Eq. JHJ without corrections for quan- 
tum effects. The dotted lines are corrected for quantum effects in 
accordance with Sec. 13.21 



f in such a way, that the classical ion crystal at F — f has 
approximately the same root-mean-square displacement as 
the quantum system at a given F and T/T p 



r = r 



1 + l T >-i 



, -1/2 



4 T 2 ', 



2,2-1-1/2 

1 + 0.013^ 



(10) 



Here u_i « 2.7986 and tt_2 ~ 12.973 are moments of 
the phonon spectrum (see iBaiko et al.l l200ll . and the val- 
ues given are for a body centered cubic crystal). We assume 
that r = F substituted in Eqs. ((6} and (|8j| provide qualita- 
tively correct results not only for classical crystals (T > T p ) 
but also for quantum crystals. In Sec. [4] it will be shown that 
the quantum corrections are typically not very large. 



4 DISCUSSION 

The knowledge of the parameters [Eq. P|] of the Zhurkov 
model of strength [Eq. |6|] allow us to estimate the breaking 
stress for timescales ~ 1 s to 1 year which are the most inter- 
esting for astrophysical applications. Direct MD simulations 
for such timescales are impossible because of the very small 
dynamical timescale of the ions ~ uj~ 1 ~ 10" 20 s. One would 
need at least a fewxlO 20 MD steps to simulate one second, 
but one time step takes approximately 10 core-seconds of 
computer time. 

Our estimates of the long time breaking stresses are 
shown in Fig. [2] for 56 Fe matter at a density 10 9 gem -3 . 
For each of 5 temperatures (T = 1 x 10 6 , 2 X 10 7 , 5 X 10 7 , 
1 x 10 8 and 2 X 10 s K) the figure contains two lines: the solid 
line corresponds to the Zhurkov model for classical crystals 
[Eqs. ((6| and ©] and the dotted line includes corrections 
for quantum effects as described in Sec. 13.21 

One can see a strong dependence of the breaking stress 
on the temperature. For the highest temperature T — 2 x 



10 K the estimated breaking stress almost vanishes on the 
timescale of a few years. This can limit the lifetime of moun- 
tains, supported by elasticity, on neutron stars with warm 
crust. For lower temperatures the breaking stress is signif- 
icantly larger and the T-dependence becomes weaker - the 
breaking stresses at T — 10 8 and T = 10 6 K differ by less 
than a factor of two. Finally, for T < 10 6 K the breaking 
stress is almost independent of temperature. The timescale 
dependence is strong for high temperatures T ~ 2 x 10 8 K 
and almost vanishes for low temperatures T < 2 x 10 7 K. 

In our model (Sec. 13. 2[) corrections for quantum effects 
are not very large (compare solid and dotted lines in Fig. 
0). The reason is simple. At large temperatures the ions are 
almost classical, and quantum corrections are small. To de- 
scribe low temperature limit let us estimate breaking stress 
(j b w U /N + A(T b , where correction A(j b w — ln(roj p ) / (TN) 
is small at low temperature. In our model the quantum cor- 
rections are described by decreasing of the effective coupling 
parameter T, and affect only the correction therm Aer b . 

We should note, that our estimates for 1 s to 1 year 
timescales are based on the extrapolation of the MD data 
of more than ten orders of magnitude in time. However our 
data only span three orders of magnitude in time (10 -7 < 
v/ujp < 10~ 4 ). So the extrapolation results should be taken 
with caution, especially for low coupling V < 250 where the 
dependence of the breaking stress on time is most significant. 
For large enough T > 800 the breaking stress extrapolated 
to timescales of a few years is only slightly lower than at 
timescales achieved in our MD simulations (compare Fig. [T] 
and the right axis of Fig. [2}. We expect that the extrapola- 
tion for such strong coupling is more reliable. However, we 
can not exclude that there are instabilities of the deformed 
crystal which have large enough timescales to be insignifi- 
cant in our MD simulations but could reduce the breaking 
stress for few second timescales. In any case, the breaking 
stress at long timescales can not exceed the breaking stress 
a™ 1 " 1 obtained in our longest simulations. The corresponding 
values can be fitted as 



max 
0"b 



0.0195 



1.27 

r- 7i 



(ii) 



This fit provides an upper limit for the long time breaking 
stress. 



5 CONCLUSIONS 

We have performed extensive MD studies of the breaking 
stress of neutron star crust. Our results are in good agree- 
ment with the Zhurkov model of strength, and we have de- 
termined the corresponding parameters, Eqs. ^ and ©. 
We apply this parametrization to estimate the breaking 
stress for very long timescales ~ 1 s to 1 year, which are 
too large for direct MD studies. We demonstrate that for a 
coupling parameter V < 200 neutron star crust matter can 
break under very small stress, if it is applied for a few years. 
This result is based on an extrapolation of over ten orders of 
magnitude in timescale and should be treated with caution. 
We construct a qualitative model to include the influence 
of quantum corrections on breaking stress (Sec. I3.2[) . We 
show that these corrections are small (solid and dotted lines 
in Fig. [2]). We also present a fit for the upper limit on the 
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breaking stress for long times, Eq. which is based on 

the breaking stresses for the slowest deformations directly 
obtained in our MD simulations. Our results can be used to 
estimate the lifetime of mountains on neutron stars and to 
obtain the breaking stress for the relevant timescales of other 
astrophysical phenomena associated with crust breaking. 

In this paper we concentrated on shear deformations 
of a perfect body centered cubic crystal along the lattice 
planes. Realistic neutron star crust can contain impurities, 
grain boundaries and the deformation may not be aligned 
with th e lattice. These effects we re found to be unimpor- 
tant bv lHorowitz fc Kadaul {2009), but more detailed stud- 
ies are need. In this paper we also considered only one ratio 
of screening length to ion sphere radius A c /a = 1.75, but 
this ratio depends on the composition and can affect the 
breaking stress. In addition, breaking stress for the tension 
deformations can be important. We plan to study these ef- 
fects in the near future. 
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